import numpy as np

def normal(X):
    total = sum(X)
    return [x / total for x in X]

def printf(result):
    print("[", end=' ')
    for num in result:
        print('%.3f' % num, end=' ')
    print("]")


PA=[0.3,0.3,0.2,0.2]
PB=[0.4,0.4,0.1,0.1]
PC=[0.2,0.2,0.3,0.3]

PS=[[[0.2,0.6,0.2],[0.1,0.3,0.6],[0.05,0.2,0.75],[0.01,0.1,0.89]],
[[0.6,0.3,0.1],[0.2,0.6,0.2],[0.1,0.3,0.6],[0.05,0.2,0.75]],
[[0.75,0.2,0.05],[0.6,0.3,0.1],[0.2,0.6,0.2],[0.1,0.3,0.6]],
[[0.89,0.1,0.01],[0.75,0.2,0.05],[0.6,0.3,0.1],[0.2,0.6,0.2]]]

# 精确求解方法
def direct_cal():
    res=[0,0,0]
    for XA in range(4):
        for XB in range(4):
            for XC in range(4):
                for sBC in range(3):
                    res[sBC]+=PA[XA]*PB[XB]*PC[XC]*PS[XA][XB][0]*PS[XA][XC][1]*PS[XB][XC][sBC]
    return normal(res)   #normal(X)=X/sum(X)将计数变成概率
    # return res

# 拒绝采样方法
def reject_sampling():
    n=5000
    res=[0,0,0]
    for i in range(n):
        XA=np.random.choice(4,p=PA)
        XB=np.random.choice(4,p=PB)
        XC=np.random.choice(4,p=PC)
        sAB=np.random.choice(3,p=PS[XA][XB])
        sAC=np.random.choice(3,p=PS[XA][XC])
        sBC=np.random.choice(3,p=PS[XB][XC])
        if sAB==0 and sAC==1:
            res [sBC]+=1
    return normal(res)
    # return res

# 似然加权采样方法
def likehood_weighting():
    n=5000
    res=[0,10,0]
    for i in range(n):
        w=1
        XA=np.random.choice(4,p=PA)
        XB=np.random.choice(4,p=PB)
        XC=np.random.choice(4,p=PC)

        w=w*PS[XA][XB][0] #sAB加权
        w=w*PS[XA][XC][1] #sAC加权
        sBC=np.random.choice(3,p=PS[XB][XC])
        res[sBC]+=w
    return normal(res)
    # return res


# Gibbs采样方法
def Gibs():
    n=4999
    res=[0,0,0]
    XA,XB,XC,sAB,sAC,sBC=0,0,1,0,1,1
    for k in range(n):
        _PA=normal([PA[i]*PS[i][XB][sAB]*PS[i][XC][sAC]
                   for i in range(4)])
        XA=np.random.choice(4,p=PA)

        _PB=normal([PB[i]*PS[XA][i][sAB]*PS[i][XC][sBC]
                   for i in range(4)])
        XB=np.random.choice(4,p=PB)

        _PC = normal([PC[i] * PS[XA][i][sAC] * PS[XB][i][sBC]
                     for i in range(4)])
        XC = np.random.choice(4, p=PC)

        sBC = np.random.choice(3, p=PS[XB][XC])
        res[sBC] += 1
    return normal(res)
    # return res

printf(direct_cal())
printf(reject_sampling())
printf(likehood_weighting())
printf(Gibs())


